Load the necessary libraries
library(car) #for regression diagnostics
library(broom) #for tidy output
library(ggfortify) #for model diagnostics
library(sjPlot) #for outputs
library(knitr) #for kable
library(effects) #for partial effects plots
library(ggeffects) #for partial effects plots
library(emmeans) #for estimating marginal means
library(MASS) #for glm.nb
library(MuMIn) #for AICc
library(tidyverse) #for data wrangling
library(broom.mixed)
library(nlme) #for lme
library(lme4) #for lmer
library(lmerTest) #for Satterthwaite's p-values
library(glmmTMB) #for glmmTMB
library(DHARMa) #for residuals and diagnostics
library(performance) #for diagnostic plots
library(see) #for diagnostic plots
Starlings
Sampling design
Format of starling_full.RSV data files
| SITUATION | MONTH | MASS | BIRD |
|---|---|---|---|
| tree | Nov | 78 | tree1 |
| .. | .. | .. | .. |
| nest-box | Nov | 78 | nest-box1 |
| .. | .. | .. | .. |
| inside | Nov | 79 | inside1 |
| .. | .. | .. | .. |
| other | Nov | 77 | other1 |
| .. | .. | .. | .. |
| tree | Jan | 85 | tree1 |
| .. | .. | .. | .. |
| SITUATION | Categorical listing of roosting situations (tree, nest-box, inside or other) |
| MONTH | Categorical listing of the month of sampling. |
| MASS | Mass (g) of starlings. |
| BIRD | Categorical listing of individual bird repeatedly sampled. |
This is a split-plot (or repeated measures) design. The individual birds are the blocks, the Situation is the between block effect and the Month is the within block effect. Repeated measures analyses involve a within block effect that represents time (in this case Month). Since it is not possible to randomize the order of time, repeated measures designs have the potential for the residuals to be autocorrelated. That is, rather than being independent, residuals from observations that are closer in time, tend to be more similar (correlated) than the residuals associated with observations that are further apart in time.
That said, with only two time points, autocorrelation is not possible.
starling = read_csv('../public/data/starling_full.csv', trim_ws=TRUE)
glimpse(starling)
## Rows: 80
## Columns: 5
## $ MONTH <chr> "Nov", "Nov", "Nov", "Nov", "Nov", "Nov", "Nov", "Nov", "No…
## $ SITUATION <chr> "tree", "tree", "tree", "tree", "tree", "tree", "tree", "tr…
## $ subjectnum <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 1, 2, 3, 4, 5, 6, 7, 8, 9, 1…
## $ BIRD <chr> "tree1", "tree2", "tree3", "tree4", "tree5", "tree6", "tree…
## $ MASS <dbl> 78, 88, 87, 88, 83, 82, 81, 80, 80, 89, 78, 78, 85, 81, 78,…
Lets prepare the data:
starling = starling %>% mutate(MONTH=factor(MONTH,levels=c('Nov','Jan')),
SITUATION=factor(SITUATION),
BIRD=factor(BIRD)
)
Model formula: \[ y_i \sim{} \mathcal{N}(\mu_i, \sigma^2)\\ \mu_i =\boldsymbol{\beta} \bf{X_i} + \boldsymbol{\gamma} \bf{Z_i} \]
where \(\boldsymbol{\beta}\) and \(\boldsymbol{\gamma}\) are vectors of the fixed and random effects parameters respectively and \(\bf{X}\) is the model matrix representing the overall intercept and effects of roosting situation and month on starling mass. \(\bf{Z}\) represents a cell means model matrix for the random intercepts associated with individual birds.
ggplot(starling, aes(y=MASS, x=MONTH)) +
geom_boxplot() +
facet_grid(~SITUATION)
## Better still
ggplot(starling, aes(y=MASS, x=MONTH, group=BIRD)) +
geom_point() +
geom_line() +
facet_grid(~SITUATION)
Conclusions:
##Multiplicative and additive models
starling.lme1 = lme(MASS ~ MONTH*SITUATION, random=~1|BIRD, data=starling, method='ML')
starling.lme2 <- update(starling.lme1, .~.-MONTH:SITUATION)
AICc(starling.lme1, starling.lme2)
Conclusions:
starling.lme2a <- update(starling.lme2, method='REML')
starling.lme2b <- update(starling.lme2a, random=~MONTH|BIRD)
AICc(starling.lme2a, starling.lme2b)
Conclusions:
##Multiplicative and additive models
starling.lmer1 = lmer(MASS ~ MONTH*SITUATION +(1|BIRD), data=starling, REML=TRUE)
starling.lmer2 <- update(starling.lmer1, .~.-MONTH:SITUATION)
AICc(starling.lmer1, starling.lmer2)
Conclusions:
starling.lmer2a <- update(starling.lmer2, REML=TRUE)
#starling.lmer2b <- update(starling.lmer2a, ~ . - (1|BIRD) + (MONTH|BIRD))
#AICc(starling.lmer2a, starling.lmer2b)
Conclusions:
##Multiplicative and additive models
starling.glmmTMB1 = glmmTMB(MASS ~ MONTH*SITUATION +(1|BIRD), data=starling, REML=FALSE)
starling.glmmTMB2 = glmmTMB(MASS ~ MONTH+SITUATION +(1|BIRD), data=starling, REML=FALSE)
# OR
starling.glmmTMB2 <- update(starling.glmmTMB1, .~.-MONTH:SITUATION)
AICc(starling.glmmTMB1, starling.glmmTMB2)
Conclusions:
starling.glmmTMB2a <- update(starling.glmmTMB2, REML=TRUE)
starling.glmmTMB2b = glmmTMB(MASS ~ MONTH+SITUATION +(MONTH|BIRD), data=starling, REML=TRUE)
# OR
starling.glmmTMB2b <- update(starling.glmmTMB2a, ~ . - (1|BIRD) + (MONTH|BIRD))
AICc(starling.glmmTMB2a, starling.glmmTMB2b)
Conclusions:
plot_grid(plot_model(starling.lme2a, type='diag')[-2])
Conclusions:
performance::check_model(starling.lme2a)
## Could not compute standard errors from random effects for diagnostic plot.
## Homogeneity of variance could not be computed. Cannot extract residual variance from objects of class 'lme'.
Conclusions:
#starling.resid = simulateResiduals(starling.lme2a, plot=TRUE)
plot_grid(plot_model(starling.lmer2a, type='diag')[-2])
Conclusions:
performance::check_model(starling.lmer2a)
Conclusions:
starling.resid = simulateResiduals(starling.lmer2a, plot=TRUE)
Conclusions:
plot_grid(plot_model(starling.glmmTMB2a, type='diag')[-2])
Conclusions:
performance::check_model(starling.glmmTMB2a)
Conclusions:
starling.resid = simulateResiduals(starling.glmmTMB2a, plot=TRUE)
Conclusions:
plot_model(starling.lme2a, type='eff', terms=c('SITUATION', 'MONTH'))
plot_model(starling.lme2a, type='est')
plot(allEffects(starling.lme2a) , multiline=TRUE, ci.style='bars')
ggpredict(starling.lme2a, terms=c('SITUATION', 'MONTH')) %>% plot()
#ggemmeans(starling.lme2a, ~SITUATION*MONTH) %>% plot()
plot_model(starling.lmer2a, type='eff', terms=c('SITUATION', 'MONTH'))
plot_model(starling.lmer2a, type='est')
plot(allEffects(starling.lmer2a) , multiline=TRUE, ci.style='bars')
ggpredict(starling.lmer2a, terms=c('SITUATION', 'MONTH')) %>% plot()
ggemmeans(starling.lmer2a, ~SITUATION*MONTH) %>% plot()
plot_model(starling.glmmTMB2a, type='eff', terms=c('SITUATION', 'MONTH'))
plot_model(starling.glmmTMB2a, type='est')
plot(allEffects(starling.glmmTMB2a) , multiline=TRUE, ci.style='bars')
ggpredict(starling.glmmTMB2a, terms=c('SITUATION', 'MONTH')) %>% plot()
ggemmeans(starling.glmmTMB2a, ~SITUATION*MONTH) %>% plot()
summary(starling.lme2a)
## Linear mixed-effects model fit by REML
## Data: starling
## AIC BIC logLik
## 456.198 472.4204 -221.099
##
## Random effects:
## Formula: ~1 | BIRD
## (Intercept) Residual
## StdDev: 0.7571596 4.110025
##
## Fixed effects: MASS ~ MONTH + SITUATION
## Value Std.Error DF t-value p-value
## (Intercept) 78.85 1.0550347 39 74.73688 0.0000
## MONTHJan 9.10 0.9190296 39 9.90175 0.0000
## SITUATIONnest-box 1.40 1.3430893 36 1.04237 0.3042
## SITUATIONother -3.60 1.3430893 36 -2.68039 0.0110
## SITUATIONtree 3.80 1.3430893 36 2.82930 0.0076
## Correlation:
## (Intr) MONTHJ SITUATION- SITUATIONth
## MONTHJan -0.436
## SITUATIONnest-box -0.637 0.000
## SITUATIONother -0.637 0.000 0.500
## SITUATIONtree -0.637 0.000 0.500 0.500
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -1.7295431 -0.6541972 -0.0445886 0.7190703 2.0485226
##
## Number of Observations: 80
## Number of Groups: 40
Conclusions:
fixef(starling.lme2a)
## (Intercept) MONTHJan SITUATIONnest-box SITUATIONother
## 78.85 9.10 1.40 -3.60
## SITUATIONtree
## 3.80
tidy(starling.lme2a, conf.int=TRUE)
tidy(starling.lme2a, conf.int=TRUE) %>% kable
| effect | group | term | estimate | std.error | df | statistic | p.value | conf.low | conf.high |
|---|---|---|---|---|---|---|---|---|---|
| fixed | fixed | (Intercept) | 78.8500000 | 1.0550347 | 39 | 74.736876 | 0.0000000 | 76.715991 | 80.9840092 |
| fixed | fixed | MONTHJan | 9.1000000 | 0.9190296 | 39 | 9.901749 | 0.0000000 | 7.241087 | 10.9589128 |
| fixed | fixed | SITUATIONnest-box | 1.4000000 | 1.3430893 | 36 | 1.042373 | 0.3041886 | -1.323911 | 4.1239114 |
| fixed | fixed | SITUATIONother | -3.6000000 | 1.3430893 | 36 | -2.680388 | 0.0110226 | -6.323911 | -0.8760886 |
| fixed | fixed | SITUATIONtree | 3.8000000 | 1.3430893 | 36 | 2.829298 | 0.0075798 | 1.076089 | 6.5239114 |
| ran_pars | BIRD | sd_(Intercept) | 0.7571596 | NA | NA | NA | NA | 0.004007 | 143.0736430 |
| ran_pars | Residual | sd_Observation | 4.1100253 | NA | NA | NA | NA | NA | NA |
Conclusions:
# warning this is only appropriate for html output
sjPlot::tab_model(starling.lme2a, show.se=TRUE, show.aic=TRUE)
| Â | MASS | |||
|---|---|---|---|---|
| Predictors | Estimates | std. Error | CI | p |
| (Intercept) | 78.85 | 1.06 | 76.72 – 80.98 | <0.001 |
| MONTH [Jan] | 9.10 | 0.92 | 7.24 – 10.96 | <0.001 |
| SITUATION [nest-box] | 1.40 | 1.34 | -1.32 – 4.12 | 0.304 |
| SITUATION [other] | -3.60 | 1.34 | -6.32 – -0.88 | 0.011 |
| SITUATION [tree] | 3.80 | 1.34 | 1.08 – 6.52 | 0.008 |
| Random Effects | ||||
| σ2 | 16.89 | |||
| τ00 BIRD | 0.57 | |||
| ICC | 0.03 | |||
| N BIRD | 40 | |||
| Observations | 80 | |||
| Marginal R2 / Conditional R2 | 0.618 / 0.630 | |||
| AIC | 456.198 | |||
summary(starling.lmer2a)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: MASS ~ MONTH + SITUATION + (1 | BIRD)
## Data: starling
##
## REML criterion at convergence: 442.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.72954 -0.65420 -0.04459 0.71907 2.04852
##
## Random effects:
## Groups Name Variance Std.Dev.
## BIRD (Intercept) 0.5733 0.7572
## Residual 16.8923 4.1100
## Number of obs: 80, groups: BIRD, 40
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 78.850 1.055 52.189 74.737 < 2e-16 ***
## MONTHJan 9.100 0.919 39.000 9.902 3.38e-12 ***
## SITUATIONnest-box 1.400 1.343 36.000 1.042 0.30419
## SITUATIONother -3.600 1.343 36.000 -2.680 0.01102 *
## SITUATIONtree 3.800 1.343 36.000 2.829 0.00758 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) MONTHJ SITUATION- SITUATIONth
## MONTHJan -0.436
## SITUATIONn- -0.637 0.000
## SITUATIONth -0.637 0.000 0.500
## SITUATIONtr -0.637 0.000 0.500 0.500
Conclusions:
confint(starling.lmer2a)
## 2.5 % 97.5 %
## .sig01 0.000000 2.353286
## .sigma 3.491034 4.763279
## (Intercept) 76.844539 80.855461
## MONTHJan 7.306261 10.893739
## SITUATIONnest-box -1.136730 3.936730
## SITUATIONother -6.136730 -1.063270
## SITUATIONtree 1.263270 6.336730
tidy(starling.lmer2a, conf.int=TRUE)
tidy(starling.lmer2a, conf.int=TRUE) %>% kable
| effect | group | term | estimate | std.error | statistic | df | p.value | conf.low | conf.high |
|---|---|---|---|---|---|---|---|---|---|
| fixed | NA | (Intercept) | 78.8500000 | 1.0550347 | 74.736877 | 52.18865 | 0.0000000 | 76.782170 | 80.9178301 |
| fixed | NA | MONTHJan | 9.1000000 | 0.9190296 | 9.901749 | 38.99999 | 0.0000000 | 7.298735 | 10.9012649 |
| fixed | NA | SITUATIONnest-box | 1.4000000 | 1.3430893 | 1.042373 | 35.99999 | 0.3041885 | -1.232407 | 4.0324066 |
| fixed | NA | SITUATIONother | -3.6000000 | 1.3430893 | -2.680388 | 35.99999 | 0.0110226 | -6.232407 | -0.9675934 |
| fixed | NA | SITUATIONtree | 3.8000000 | 1.3430893 | 2.829298 | 35.99999 | 0.0075798 | 1.167593 | 6.4324066 |
| ran_pars | BIRD | sd__(Intercept) | 0.7571592 | NA | NA | NA | NA | NA | NA |
| ran_pars | Residual | sd__Observation | 4.1100253 | NA | NA | NA | NA | NA | NA |
Conclusions:
# warning this is only appropriate for html output
sjPlot::tab_model(starling.lmer2a, show.se=TRUE, show.aic=TRUE)
| Â | MASS | |||
|---|---|---|---|---|
| Predictors | Estimates | std. Error | CI | p |
| (Intercept) | 78.85 | 1.06 | 76.78 – 80.92 | <0.001 |
| MONTH [Jan] | 9.10 | 0.92 | 7.30 – 10.90 | <0.001 |
| SITUATION [nest-box] | 1.40 | 1.34 | -1.23 – 4.03 | 0.297 |
| SITUATION [other] | -3.60 | 1.34 | -6.23 – -0.97 | 0.007 |
| SITUATION [tree] | 3.80 | 1.34 | 1.17 – 6.43 | 0.005 |
| Random Effects | ||||
| σ2 | 16.89 | |||
| τ00 BIRD | 0.57 | |||
| ICC | 0.03 | |||
| N BIRD | 40 | |||
| Observations | 80 | |||
| Marginal R2 / Conditional R2 | 0.618 / 0.630 | |||
| AIC | 456.198 | |||
summary(starling.glmmTMB2a)
## Family: gaussian ( identity )
## Formula: MASS ~ MONTH + SITUATION + (1 | BIRD)
## Data: starling
##
## AIC BIC logLik deviance df.resid
## 456.2 472.9 -221.1 442.2 78
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## BIRD (Intercept) 0.5733 0.7572
## Residual 16.8923 4.1100
## Number of obs: 80, groups: BIRD, 40
##
## Dispersion estimate for gaussian family (sigma^2): 16.9
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 78.850 1.055 74.74 < 2e-16 ***
## MONTHJan 9.100 0.919 9.90 < 2e-16 ***
## SITUATIONnest-box 1.400 1.343 1.04 0.29724
## SITUATIONother -3.600 1.343 -2.68 0.00735 **
## SITUATIONtree 3.800 1.343 2.83 0.00467 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Conclusions:
confint(starling.glmmTMB2a)
## 2.5 % 97.5 % Estimate
## cond.(Intercept) 76.782170256 80.9178297 78.8500000
## cond.MONTHJan 7.298735065 10.9012649 9.1000000
## cond.SITUATIONnest-box -1.232406127 4.0324061 1.4000000
## cond.SITUATIONother -6.232406127 -0.9675939 -3.6000000
## cond.SITUATIONtree 1.167593873 6.4324061 3.8000000
## BIRD.cond.Std.Dev.(Intercept) 0.005704079 100.5046556 0.7571568
tidy(starling.glmmTMB2a, conf.int=TRUE)
tidy(starling.glmmTMB2a, conf.int=TRUE) %>% kable
| effect | component | group | term | estimate | std.error | statistic | p.value | conf.low | conf.high |
|---|---|---|---|---|---|---|---|---|---|
| fixed | cond | NA | (Intercept) | 78.8500000 | 1.0550346 | 74.736888 | 0.0000000 | 76.7821703 | 80.9178297 |
| fixed | cond | NA | MONTHJan | 9.1000000 | 0.9190296 | 9.901748 | 0.0000000 | 7.2987351 | 10.9012649 |
| fixed | cond | NA | SITUATIONnest-box | 1.4000000 | 1.3430890 | 1.042373 | 0.2972387 | -1.2324061 | 4.0324061 |
| fixed | cond | NA | SITUATIONother | -3.6000000 | 1.3430890 | -2.680388 | 0.0073537 | -6.2324061 | -0.9675939 |
| fixed | cond | NA | SITUATIONtree | 3.8000000 | 1.3430890 | 2.829299 | 0.0046650 | 1.1675939 | 6.4324061 |
| ran_pars | cond | BIRD | sd__(Intercept) | 0.7571568 | NA | NA | NA | 0.0057041 | 100.5046556 |
| ran_pars | cond | Residual | sd__Observation | 4.1100254 | NA | NA | NA | 0.0057041 | 100.5046556 |
Conclusions:
# warning this is only appropriate for html output
sjPlot::tab_model(starling.glmmTMB2a, show.se=TRUE, show.aic=TRUE)
| Â | MASS | |||
|---|---|---|---|---|
| Predictors | Estimates | std. Error | CI | p |
| (Intercept) | 78.85 | 1.06 | 76.78 – 80.92 | <0.001 |
| MONTH [Jan] | 9.10 | 0.92 | 7.30 – 10.90 | <0.001 |
| SITUATION [nest-box] | 1.40 | 1.34 | -1.23 – 4.03 | 0.297 |
| SITUATION [other] | -3.60 | 1.34 | -6.23 – -0.97 | 0.007 |
| SITUATION [tree] | 3.80 | 1.34 | 1.17 – 6.43 | 0.005 |
| Random Effects | ||||
| σ2 | 16.89 | |||
| τ00 BIRD | 0.57 | |||
| ICC | 0.03 | |||
| N BIRD | 40 | |||
| Observations | 80 | |||
| Marginal R2 / Conditional R2 | 0.618 / 0.630 | |||
| AIC | 456.198 | |||
emmeans(starling.lme2a, pairwise~SITUATION)
## $emmeans
## SITUATION emmean SE df lower.CL upper.CL
## inside 83.4 0.95 39 81.5 85.3
## nest-box 84.8 0.95 36 82.9 86.7
## other 79.8 0.95 36 77.9 81.7
## tree 87.2 0.95 36 85.3 89.1
##
## Results are averaged over the levels of: MONTH
## Degrees-of-freedom method: containment
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## inside - (nest-box) -1.4 1.34 36 -1.042 0.7259
## inside - other 3.6 1.34 36 2.680 0.0515
## inside - tree -3.8 1.34 36 -2.829 0.0364
## (nest-box) - other 5.0 1.34 36 3.723 0.0036
## (nest-box) - tree -2.4 1.34 36 -1.787 0.2961
## other - tree -7.4 1.34 36 -5.510 <.0001
##
## Results are averaged over the levels of: MONTH
## Degrees-of-freedom method: containment
## P value adjustment: tukey method for comparing a family of 4 estimates
r.squaredGLMM(starling.lme2a)
## R2m R2c
## [1,] 0.6178293 0.6303737
## # R2 for Mixed Models
##
## Conditional R2: 0.630
## Marginal R2: 0.618
emmeans(starling.lmer2a, pairwise~SITUATION)
## $emmeans
## SITUATION emmean SE df lower.CL upper.CL
## inside 83.4 0.95 36 81.5 85.3
## nest-box 84.8 0.95 36 82.9 86.7
## other 79.8 0.95 36 77.9 81.7
## tree 87.2 0.95 36 85.3 89.1
##
## Results are averaged over the levels of: MONTH
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## inside - (nest-box) -1.4 1.34 36 -1.042 0.7259
## inside - other 3.6 1.34 36 2.680 0.0515
## inside - tree -3.8 1.34 36 -2.829 0.0364
## (nest-box) - other 5.0 1.34 36 3.723 0.0036
## (nest-box) - tree -2.4 1.34 36 -1.787 0.2961
## other - tree -7.4 1.34 36 -5.510 <.0001
##
## Results are averaged over the levels of: MONTH
## Degrees-of-freedom method: kenward-roger
## P value adjustment: tukey method for comparing a family of 4 estimates
r.squaredGLMM(starling.lmer2a)
## R2m R2c
## [1,] 0.6178293 0.6303737
## # R2 for Mixed Models
##
## Conditional R2: 0.630
## Marginal R2: 0.618
emmeans(starling.glmmTMB2a, pairwise~SITUATION)
## $emmeans
## SITUATION emmean SE df lower.CL upper.CL
## inside 83.4 0.95 78 81.5 85.3
## nest-box 84.8 0.95 78 82.9 86.7
## other 79.8 0.95 78 77.9 81.7
## tree 87.2 0.95 78 85.3 89.1
##
## Results are averaged over the levels of: MONTH
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## inside - (nest-box) -1.4 1.34 78 -1.042 0.7252
## inside - other 3.6 1.34 78 2.680 0.0436
## inside - tree -3.8 1.34 78 -2.829 0.0296
## (nest-box) - other 5.0 1.34 78 3.723 0.0021
## (nest-box) - tree -2.4 1.34 78 -1.787 0.2873
## other - tree -7.4 1.34 78 -5.510 <.0001
##
## Results are averaged over the levels of: MONTH
## P value adjustment: tukey method for comparing a family of 4 estimates
r.squaredGLMM(starling.glmmTMB2a)
## R2m R2c
## [1,] 0.6178294 0.6303736
## # R2 for Mixed Models
##
## Conditional R2: 0.630
## Marginal R2: 0.618
newdata = emmeans(starling.lme2a, ~SITUATION*MONTH) %>%
as.data.frame
ggplot(newdata, aes(y=emmean, x=SITUATION, fill=MONTH)) +
geom_pointrange(aes(ymin=lower.CL, ymax=upper.CL), shape=21,
position=position_dodge(0.2)) +
theme_classic()
newdata = emmeans(starling.lmer2a, ~SITUATION*MONTH) %>%
as.data.frame
ggplot(newdata, aes(y=emmean, x=SITUATION, fill=MONTH)) +
geom_pointrange(aes(ymin=lower.CL, ymax=upper.CL), shape=21,
position=position_dodge(0.2)) +
theme_classic()
newdata = emmeans(starling.glmmTMB2a, ~SITUATION*MONTH) %>%
as.data.frame
ggplot(newdata, aes(y=emmean, x=SITUATION, fill=MONTH)) +
geom_pointrange(aes(ymin=lower.CL, ymax=upper.CL), shape=21,
position=position_dodge(0.2)) +
theme_classic()